%% Matlab codes to generate the main results in "Asymmetric Network Connectedness of Fears" by Barunik, Bevilacqua, Tunaru.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 1) IV Data and plot.

%Plot IV Series -- CitiGroup 
%IV_Data.mat contains the main dataset of the paper: the total and decomposed VIX series and the dates column. 
load('IV_Data.mat')
Date = datenum(Date);

h = figure;
plot(Date,IV_Data{:,'CITIVIX'},'DisplayName','CITIVIX','Color',[0 0 0]);hold on;
plot(Date,IV_Data{:,'CITIVIXM'},'DisplayName','CITIVIXM','LineStyle','-.',...
    'Color',[0.8 0.8 0.8]);hold on;
plot(Date,IV_Data{:,'CITIVIXP'},'DisplayName','CITIVIXP','LineStyle','--',...
    'Color',[0.501960784313725 0.501960784313725 0.501960784313725]);hold off;
recessionplot
%title('CITI IVs')
datetick('x')
grid on
hold off;
set(gca,'LooseInset',get(gca,'TightInset'));
set(gcf,'position',[2,2,500,200])
set(h,'Units','Inches');
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h,'Figure1.pdf','-dpdf')

clc
clear all;
%% 2) DY Network Index -- Implied Volatility Asymmetry US. 

% This second part of the code is divided into three parts and the codes
% compute the network connectedness with respect to the total IV, negative IV, 
% and positive IV, all in both static framework and in dynamics.

% we split the IV_Data.mat whole data set in the different IV classes and
% converted as double for ease to work with a single group of IVs. 

% for this replication code, we load all of them, and show how the code works 
% with respect to the IV total.  

%load .mat data files.
clc
clear all
%IV TOT
load('IVTOT_Fin_Sect.mat'); 
%IV minus
load('IVMIN_Fin_Sect.mat');
%IV plus
load('IVPLUS_Fin_Sect.mat');

%load dates
load('Date_IV');
%% a) Compute the Total Connectedness Network of Volatility

%Data to be input in the network here -- e.g. input total VIX series -- "IVTOT".
Dat_ = IVTOT;  %replace IV series to study here
%convert them in log volatilities.
Dat_=log(Dat_);
%Date transformation
DateNum = datenum(Date);
% Inputs to the SpillOver_DY function are: SpillOver_DY(VAR_lag, nMA_lags, data, plot_)
%parameters set as discussed in the paper. 
nObs = length(Dat_);
[ Out_, Spill_indx, IRF_g, V_decom_g ] = DYFunctReplica(4, 12, Dat_, 0);
%VD will give us the variance decomposition matrix, Out is the FROM and TO
%statistics, Spill_indx will give us the total spillover index value.
VD=V_decom_g*100;
Out=Out_*100;
O=Out(:,2)';

% this will reproduce Table 1. replace data in line 58 to construct the 
% connectedness variance decomposition table with respect to the decomposed
%volatility measures (Table 2).

%% b) Rolling Window calculations,
%
WL      = 200;                      %Set window length, in our case = 200 days. 
x_start = [1:1:nObs-WL+1 ]';
x_slut  = [WL:1:nObs ]';
n_iter  = length(x_start);
indx_   = zeros(n_iter,1);
Date(1:WL-1) = [ ];
DateNum(1:WL-1) = [ ];

%% Total Spillover Index
h = waitbar(0,'Iterating, please wait...');
for ( j=1:n_iter )
    waitbar(j/n_iter,h)
    [ Out_, Spill_indx, IRF_g, V_decom_g ] = DYFunctReplica(4, 12, Dat_(x_start(j,1):x_slut(j,1),:), 0); 
    indx_(j,1) = Spill_indx(1,1);
end
close(h)

%this will give us the Total Fear Connectedness Index over time (dates - window lenght)
%replicate Figure 2 of the paper and save it. 
indx = indx_*100;

h = figure
plot(Date, indx)
set(gca,'LooseInset',get(gca,'TightInset'));
set(gcf,'position',[2,2,500,200])
set(h,'Units','Inches');
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h,'Figure2.pdf','-dpdf')

%(we saved it already "C_IV_TOT_200.mat", and dates "C_IV_Date.mat".
% uncomment lines 112, 113 to overwrite.

%filename = 'C_IV_TOT_200.mat';
%save(filename, 'indx')

%replace data in line 58 to construct the connectedness measure for
%the decomposed volatility measures and reproduce Figure 3. 
%We saved them already as:
% 1) 'C_IV_MIN_200.mat'
% 2) C_IV_PLUS_200.mat'
% and work with them in the next section. 

%% c) AFC

%Compute AFC - Asymmetric Fear Connectedness measure - As the difference between
%%C_IV+ - C_IV-. When AFC is positive implies that the connectedness coming from
%%IV+ are larger than the ones coming from IV-, the opposite is true when
%%AFC is negative. See Section 3 of the paper for more details. 

%We saved the output of the previous section's code when it is computed for
% all IV series, namely total, minus and plus volatility series. We load
% them here together with the date series minus the rolling window length.
% We compute the AFC series as the difference between  C^+ and C^-. 
load('C_IV_TOT_200.mat');
load('C_IV_MIN_200.mat');
load('C_IV_PLUS_200.mat');
load('C_IV_Date.mat');
Date = datenum(Date);
AFC= C_IV_PLUS-C_IV_MIN;

% We replicate here Figure 3 in which we plot the two decomposed positive and negative
% IV connectedness measures together with their AFC in a multiplot.
h = figure 
subplot(2,1,1)
plot(Date,C_IV_TOT,'DisplayName','C VIX TOT');hold on;plot(Date,C_IV_MIN,'DisplayName','C VIX MIN');hold on;plot(Date,C_IV_PLUS,'DisplayName','C VIX_PLUS');recessionplot;hold off
title('VIX Connectedness Indexes')
datetick('x')
grid on
subplot(2,1,2)
plot(Date, AFC,'DisplayName','US Financial Sector AFC');recessionplot;
title('US Financial Sector AFC')
datetick('x')
grid on
hold off;
% Add lines
h1 = line([733660 733660],[-3.5 7.3]);
h2 = line([736876 736876],[-3.5 7.3]);
% Set properties of lines
set([h1 h2],'Color','w','LineWidth',1)
% Add a patch
%patch([733660 736876 736876 733660],[-3.5 -3.5 7.3 7.3],'c')
% The order of the "children" of the plot determines which one appears on top.
% I need to flip it here.
set(gca,'children',flipud(get(gca,'children')))
hold off
set(gca,'LooseInset',get(gca,'TightInset'));
set(gcf,'position',[2,2,1000,800])
set(h,'Units','Inches');
pos = get(h,'Position');
set(h,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
print(h,'Figure3.pdf','-dpdf')

%% d) NET IV Connectedness Measures

% We here report the code to compute the NET IV connectedness measures as
% difference between TO and FROM connectedness as detailed in the paper in
% Section 3. We also plot them for each bank in the online appendix.

%% d.1) Rolling Window calculation also for the "Directional TO others"
TO_=[ ];
TO_   = zeros(n_iter,10); 
h = waitbar(0,'Iterating, please wait...');
for ( j=1:n_iter )
    waitbar(j/n_iter,h)
    [ Out_, Spill_indx, IRF_g, V_decom_g ] = DYFunctReplica(4, 12, Dat_(x_start(j,1):x_slut(j,1),:), 0); 
    TO_(j,1) = Out_(1,2); 
    TO_(j,2) = Out_(2,2);
    TO_(j,3) = Out_(3,2); 
    TO_(j,4) = Out_(4,2);
    TO_(j,5) = Out_(5,2); 
    TO_(j,6) = Out_(6,2);
    TO_(j,7) = Out_(7,2); 
    TO_(j,8) = Out_(8,2);
     TO_(j,9) = Out_(9,2); 
    TO_(j,10) = Out_(10,2);
end
close(h)
%% d.2) Rolling Window calculation also for the "Directional FROM others"

h = waitbar(0,'Iterating, please wait...');
for ( j=1:n_iter )
    waitbar(j/n_iter,h)
    [ Out_, Spill_indx, IRF_g, V_decom_g ] = DYFunctReplica(4, 12, Dat_(x_start(j,1):x_slut(j,1),:), 0); 
    FROM_(j,1) = Out_(1,1); 
    FROM_(j,2) = Out_(2,1);
    FROM_(j,3) = Out_(3,1); 
    FROM_(j,4) = Out_(4,1);
    FROM_(j,5) = Out_(5,1); 
    FROM_(j,6) = Out_(6,1);
    FROM_(j,7) = Out_(7,1); 
    FROM_(j,8) = Out_(8,1);
    FROM_(j,9) = Out_(9,1); 
    FROM_(j,10) = Out_(10,1);
end
close(h)

%% d.3) Net Directional Connectedness 

% After computing the directional TO and FROM connectedness by running d.1 and d.2, we take the difference
% NET = TO - FROM, and this can be done for each series of VIXs, namely Total, Positive and Negative
% we have also saved them for ease (for the total VIX only). to load them just uncomment the following 4 lines.

%load('VIX_FROM_200.mat')
%load('VIX_TO_200.mat')
%load('C_IV_Date.mat')
%DateNum = datenum(Date);
%This is the formula to obtain the NET VIX connectedness from bank i to all the other banks j.
NET_ = [ ];
NET_(:,1) = TO_(:,1) - FROM_(:,1);
NET_(:,2) = TO_(:,2) - FROM_(:,2);
NET_(:,3) = TO_(:,3) - FROM_(:,3);
NET_(:,4) = TO_(:,4) - FROM_(:,4);
NET_(:,5) = TO_(:,5) - FROM_(:,5);
NET_(:,6) = TO_(:,6) - FROM_(:,6);
NET_(:,7) = TO_(:,7) - FROM_(:,7);
NET_(:,8) = TO_(:,8) - FROM_(:,8); 
NET_(:,9) = TO_(:,9) - FROM_(:,9);
NET_(:,10) = TO_(:,10) - FROM_(:,10);

% the following generate the NET connectedness plot in online appendix of the paper.
figure
subplot(5,2,1)
bar(DateNum, NET_(:,1))
title('JPM - NET')
datetick('x')
grid on
subplot(5,2,2)
bar(DateNum, NET_(:,2))
title('BAC - NET')
datetick('x')
grid on
subplot(5,2,3)
bar(DateNum, NET_(:,3))
title('WFC - NET')
datetick('x')
grid on
subplot(5,2,4)
bar(DateNum, NET_(:,4))
title('CITI - NET')
datetick('x')
grid on
subplot(5,2,5)
bar(DateNum, NET_(:,5))
title('GS - NET')
datetick('x')
grid on
subplot(5,2,6)
bar(DateNum, NET_(:,6))
title('MS - NET')
datetick('x')
grid on
subplot(5,2,7)
bar(DateNum, NET_(:,7))
title('USB - NET')
datetick('x')
grid on
subplot(5,2,8)
bar(DateNum, NET_(:,8))
title('AXP - NET')
datetick('x')
grid on
subplot(5,2,9)
bar(DateNum, NET_(:,9))
title('PNC - NET')
datetick('x')
grid on
subplot(5,2,10)
bar(DateNum, NET_(:,10))
title('BK - NET')
datetick('x')
grid on
hold off


%% e) Case Study: Goldman Sachs
% We here replicate the case study in section 5 of our paper with respect to Goldman Sachs (GS). We plot GS NET directional positive and
% NET directinal negative statistics as well as the AFC measure taken as difference between the two. We have already saved (for ease) the NET measures 
%resulting from the code in section d) when IV minus ('IVMIN_Fin_Sect.mat') and IV plus ('IVPLUS_Fin_Sect.mat') are input in line 58. 
load('NET_PLUS_200.mat');
load('NET_MIN_200.mat');
load('C_IV_Date.mat');
DateNum = datenum(Date);
AFC_NET = NET_plus - NET_minus;
figure
bar(DateNum, AFC_NET(:,5)); set(gca,'children',flipud(get(gca,'children')));hold on; plot(DateNum,NET_plus(:,5),'DisplayName','NET^+');hold on;plot(DateNum,NET_minus(:,5),'DisplayName','NET^-');recessionplot;hold off;
title('GS - NET Directional and SAM')
datetick('x')
grid on
hold off


%% 3) Out-of-sample Predictability.

%OOS Network Predictability, sample of results.
clc
clear all
input_file='Predict_Data_M';
input_sheet='C_Predict_Data';

%import date as data vector
C = xlsread(input_file, input_sheet, 'b2:b215');
Cplus = xlsread(input_file, input_sheet, 'c2:c215');
Cmin = xlsread(input_file, input_sheet, 'd2:d215');
C_ALL = [C, Cmin, Cplus];
%macro indicator  
macro = xlsread(input_file, input_sheet, 'e2:e215');

% Preliminaries
T=size(macro,1);
N=size(C_ALL,2);
R=(2007-2000)*12+1; % change in-sample period
P=T-R; % forecast evaluation period
actual=nan(P,1);
FC_HA=nan(P,1);
FC_ALL=nan(P,N+2);
FC_ALL=nan(P,2);
FC_ALL=nan(P,2);

% Computing recursive forecasts
for t=1:P;
    actual(t)=macro(R+t);
    macro_t=macro(1:R+(t-1));
    
    % change benchmark here
    %FC_HA(t)=mean(macro(1:R+(t-1)));
    %FC_HA(t) = macro(R+t-1);
    %mean over last quarter
    %FC_HA(t)=mean(macro((R+t-4):(R+t)));
    
    %autoregressive structure
    Mdl = arima(1,1,0);
    EstMdl = estimate(Mdl, macro(1:R+(t-1)));
    [FC_HA(t), FCMSE] = forecast(EstMdl,1,'Y0', macro(1:R+(t-1))); 
    C_ALL_t=C_ALL(1:R+(t-1),:);
     
    % change predictive horizons here as explained in the main text, where H = 1 + h; as well as C [1], 
    %C^+ & C^- [2:3] in lines 361,362; 
    h = 6;
    H = 1 + h;
 for i=1:N;  
        results_ALL_i_t=ols(macro_t(H:end),...
         [ones(R+(t)-H,1) C_ALL_t(1:end-h,[2:3]) macro_t(h:end-1)]);
         FC_ALL(t,i)=[1 C_ALL_t(end,[2:3]) macro_t(end)]*results_ALL_i_t.beta;
      end;
end; 


% plot of out of sample estimation vs benchmark vs actual here.
%figure; plot(actual); hold on; plot(FC_HA); hold on; plot(FC_ALL(:,1)); hold off;

% Evaluating forecasts
disp('Evaluating forecasts');
results_ALL=nan(size(FC_ALL,2),2);

for i=1:size(FC_ALL,2);
    [MSFE_adjusted_ALL_i,p_value_ALL_i]=Perform_CW_test(actual,...
        FC_HA,FC_ALL(:,i));
    results_ALL(i,:)=[MSFE_adjusted_ALL_i ...
        p_value_ALL_i];
end;
